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The Chebyshev matrix collocation method is applied to obtain the spatial modes of 
the Orr-Sommerfeld equation for Poiseuille flow and the Blausius boundary layer. The 
problem is linearized by the companion matrix technique for semi-infinite domain using 
a mapping transformation. The method can be easily adapted to problems with different 
boundary conditions requiring different transformations. 


1. INTRODUCTION 

The linear stability of hydrodynamic problems is governed by the Orr-Sommerfeld 
equation. For parallel shear flows, along with homogeneous boundary conditions, this 
equation constitutes an eigenvalue problem. For temporally evolving flows, the eigenvalue 
is the complex frequency and the problem becomes linear in the eigenvalue. Although 
for small amplitude disturbances temporal evolution is a good approximation to the 



laboratory flow, the correct physics of the problem can be obtained only through the 
calculation of spatially evolving modes. In the Orr-Sommerfeld equation, these modes 
are given in terms of the wave number and for these situations, where the complex wave 
number is the eigenvalue, the Orr-Sommerfeld equation becomes a nonlinear eigenvalue 
problem. 

The numerical solution of the Orr-Sommerfeld equation by matrix methods requires 
an efficient means of discretization. Since shear flows have steep gradients near the bound- 
aries, it is natural to seek for a discretization with clustering or streching in these regions. 
Orszag 1 ’s pioneering work demonstrated the usefulness and accuracy of Chebyshev spec- 
tral methods for solving the linear eigenvalue problems. Chebyshev methods naturally 
cluster collocation points in the vicinity of the boundaries and are more accurate than 
fourth-order finite differences on a stretched mesh with the same number of grid points. 

Several techniques to solve the temporal (linear) eigenvalue problem are discussed 
by Canuto et al . 2 , where the Chebyshev-tau method is favored over the Chebshev collo- 
cation matrix method because of the difficulties in the latter concerning the imposition of 
the boundary conditions. However, the collocation matrix method is easier to formulate 
and the boundary conditions do not pose a problem provided that they are applied near 
the boundaries, i. e. only the first and last few rows of the matrix are modified. Further- 
more, the Chebyshev-tau method requires major modifications for each new mean velocity 
profile or for a new coordinate transformation 3 . Under similar conditions, the Cheby- 
shev collocation matrix method necessitates no significant changes. Recently, Laurien* 
implemented the Chebyshev matrix collocation method to solve the temporal general- 
ized eigenvalue problem for the Blausius boundary layer and obtained very good results. 
Laurien 4 used an even number of collocation points obviating the need for the explicit 
imposition of the homogeneous boundary conditions at infinity. Note that Chebyshev 
polynomials automatically satisfy these conditions at the mid-point of the Chebyshev 
collocation domain which corresponds to infinity in the physical space. Using his tech- 
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nique with an odd number of collocation points, Biringen and Danabasoglu 5 included all 
the boundary conditions and still obtained results in excellent agreement with literature. 

The solution of nonlinear eigenvalue problems can be obtained efficiently by shooting 
methods in which a good initial guess is vital for convergence. For problems in which a 
good initial guess is not available, the matrix method presents an attractive alternative. 
Bridges and Morris 6 developed a companion matrix method to linearize the spatial eigen- 
value problem to study the stability of both channel and boundary layer flows. Their 
method incorporates the Chebyshev-tau method to discretize the governing equation. 
Once the problem is linearized by this method, the entire spectrum can be obtained via 
the QZ algorithm. Recently, Khorrami et a i. 3 used the Chebyshev collocation matrix 
method to study the temporal and spatial stability of swirling flows in enclosed domains. 
They employed the companion matrix method of Bridges and Morris 6 , however they did 
not present the eigenvector distributions corresponding to the least damped eigenvalue. 

In this work, the Chebyshev collocation matrix method is combined with the com- 
panion matrix method to solve the non-linear spatial eigenvalue problem for channel and 
Blasius boundary layer flows. The semi-infinite domain of the Blausius flow requires a 
mapping transformation and presents a challenging problem for the method under consid- 
eration. In addition to the least damped eigenvalues, we also present the corresponding 
eigenfunction distributions. 

2. GOVERNING EQUATIONS 

The stability of the parallel shear flows is governed by the Orr-Sommmerfeld equa- 
tion, 



with 

<j> = (j)' = 0 at y = ±1, (la) 

<f> = <f>' = 0 at y = 0, y — > oo. (lb) 
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Equations (la) and (lb) constitude the boundary conditions for channel and boundary 
layer flows, respectively. For spatial modes, a is the complex wave number and w is the 
real, radian frequency. Also U and U represent the base flow velocity profile and its 
second derivative, respectively. Re is the Reynolds number defined as Re = Uo8\jv for 
the boundary layer and Re = U\hjv for channel flow. Here, v is the kinematic viscosity 
of the fluid, h is the half channel height and 8\ is the displacement thickness. Velocity 
components are non-dimensionalized by the free stream velocity (Uq) and the maximum 
channel velocity (f/i) for the channel and boundary layer cases, respectively. 

Note that in the derivation of equation (1), the two-dimensional normal mode ex- 
pansion is applied using the stream function (^) formulation. Accordingly, 

V> = ^(y)e l( “ I_a,l) . (2) 

This expansion gives the following equations for the horizontal and vertical perturbation 
velocities, respectively, 

u = <f>\y)e* a *- ut \ (3) 

v = -ia<l>(y)e% a *- b,t \ (4) 

3. SOLUTION PROCEDURE 

The applied numerical procedure uses the Chebyshev-Gauss-Lobatto points for the 
y-direction discretization 


Vj = cos(nj/N), j = 0,1, ...N , (5) 

where N is the number of intervals in the domain and rf represents the coordinates of the 
collocation points in the Chebyshev domain. Recall that the Chebyshev polynomials of 
order k are defined on the interval (-1,1) by 

Tkiv) = cos(kcos~ 1 Tj), (6) 
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and 4>{rj) is expanded in Chebyshev polynomials 


N 

(T) 

/c = 0 

where a^s are the expansion coefficients. 

The solution of equation (12) depends on the basic state velocity profiles. For channel 
flow, the velocity profiles are readily given as 

U(y) = 1 — y 2 , V = 0. (8) 


Note that channel boundaries are located at y = ±1 and they coincide with the domain 
of the Chebyshev polynomials, i. e. y = rj. 

The calculation of the basic state velocity profile for the boundary layer is more 
involved and requires the solution of the Blausius boundary layer equation, 

f" + ^ = 0, (9) 

with 

m = o, 

/'( 0) = 0, (9a) 

/”( 0) = 0.332057339. 

Then, 

U(y) = f’(y)- (10) 


Note that the free stream velocity is assumed to be unity. The solution of equation (9) 
is obtained by a point-by-point Runge-Kutta-Verner method. 

The boundary layer calculations require the mapping of the physical domain onto 
the Chebyshev domain. In this work, we employ an exponential transformation given by 
Laurien 3 which maps the domain from [0, oo] to [1,0], 


V = 


e~y/ Y . 


( 11 ) 
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Here, Y (F = 20) is a mesh clustering parameter, and the boundary points are included 
in the domain satisfying the boundary conditions exactly. 

It is now convenient to write equation (1) in a different form, 

2 


D 2 — a 2 I ) - iRe< {aUI - u>I) I D - al - of/ I 


4 > = 0 , 


(12) 


where D represents the Chebyshev collocation matrix ■which gives the first derivative in 
the Chebyshev domain. The elements of D can be written as 1 

c* (-!)*+’' 


D fcj = 

Cj V k ~ Vj 

p — — Hi 


( fc T j)i 


Doo = 


2(1 -V 2 j) 
2N 2 4- 1 


= -D 


NNi 


and 


j, k = 0,1,. ..N. 

The higher order derivatives are simply obtained as powers of D, i. e. 


(13) 


Dy = D P , 


( 14 ) 


where p is the order of the derivative. In equation (12), I is the identity matrix. 

Note that equation (12) is a polynomial in a with matrix coefficients and has the 
following explicit form, 

R-4 (ct) — C 4 a 4 + C 3 T C 2 G 2 + C^a + Co, (lb) 


with 

C 4 =I, 

C 3 = iReUI, 

C 2 = -(iu>Rel + 2D 2 ), (16) 

C : = iReU" I - iReUD 2 , 

C 0 = D 4 +iReu> D 2 . 
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In the above set of equations, all the bold faced letters represent ( N + 1 ) x (JV + l) 
matrices; the last four rows of these matrices are modified for the boundary conditions. 
Note that the boundary conditions are independent of the wave number, and therefore, 
are imposed only in Co and the corresponding rows of the remaining matrices are set 
to zero. This implemantation creates a singular C 4 matrix resulting in infinite eigen- 
values. In order to remedy the infinite eigenvalue problem and to increase the accuracy 
by decreasing the number of arithmetic operations, the order of the matrices must be 
reduced. Since the zeros in the last four rows of the matrices (except Co) are to be 
preserved, the order reduction is done by simple column operations. This results in a 
4x4 upper triangular submatrix in the last four rows of the Co matrix (Fig. 1 ). If the 
boundary conditions are linearly independent, the order of the matrices is clearly reduced 
to ( N — 3) x ( N — 3) when four boundary conditions are eliminated. Here, note that each 
column switching necessitates the switching of the corresponding rows of the calculated 
eigenvector. Therefore, the original indices of the switched columns should be retained 
to decode the rows of the eigenvector once it is calculated. 

Following Bridges and Morris 6 , and recalling that the eigenvalues of the companion 
matrix are the roots of the corresponding polynomial equation, a companion matrix for 
equation (15) can be written as 

a 3 (f> N 
a 2 cj) 


{ 


/ 

-c 8 

I 

-C 2 

0 

-C a 

0 

-C 0 \ 

0 1 


/ C 4 

0 

0 

1 

0 

0 

°\ 

0 


0 

1 

0 

0 

— a 

0 

0 

1 

0 

V 

0 

0 

1 

0 ) 


V 0 

0 

0 

1/ 


> =0. 


(17) 


This equation represents a complex generalized eigenvalue problem and can be solved by 
the QZ algoritm. Note that the order of the above system is four times larger than the 
original reduced problem. According to equation (17), the physical eigenfunctions can be 
directly obtained from the last quarter of the companion matrix equation (equation (16)). 
The numerical calculations were performed on the VAX/ VMS 8550 at the University of 
Colorado at Boulder and the CYBER 205 at NASA Langley Research Center. This was 
accomplished by the use of three IMSL subroutines, GVLCG, GVCCG, CXLZ. 
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4. RESULTS, DISCUSSION AND CONCLUSION 


In this section, a comparison of our results to those of Bridges and Morris 6 and 
Jordinson 7 is presented for the stability of channel and Blausius boundary layer flows, 
respectively. 

a) Channel Flow Stability: Here, we concantrate on a test case given in Bridges 
and Morris 5 and also studied by Khorrami et a l. 3 . 

For this purpose, we set Re = 6000 and w = 0.26 which is linearly unstable for this 
flow. The first seven members of the eigenvalue spectrum are tabulated and compared 
with those of Bridges and Morris 6 in Table 1. Here increasing the order of Chebyshev 
polynomials results in noticable improvement in accuracy. However, machine and trun- 
cation errors impose a limit for the order of the polynomials that can be used: e. g. on the 
VAX/VMS 8550, increasing the order above 100 detonates the eigenvalue spectrum. This 
behaviour suggests that a trade-off exists between the order of the Chebyshev polynomials 
and the ability of the QZ algorithm to accurately solve large matrices. 

b) Blausius Boundary Layer Stability: In this section, following Jordinson 7 , we 

perform three calculations corresponding to subcritical (Re = 336, u> = 0.1297), slightly 
unstable (Re = 598, u> = 0.1201) and unstable (Re — 998, = 0.1122) cases. 

The computed least-damped eigenvalues are given in Table 2, and they are in ex- 
cellent agreement with those of Jordinson 7 . Furthermore, the eigenfunctions and their 
derivatives are plotted in Figs. (2-4). Note that increasing viscosity, i. e. decreasing Re, 
pushes the peak of the derivative distribution toward the plate. 

Another quantity of major interest is the Reynolds stress ( S ) distribution across the 
domain which is given by 7 

S = a r (4>r<t>\ - <f>r<i>i) - (18) 

where r and i refer to the real and imaginary parts of the eigenfunction or eigenvalue 
and the primes indicate differentiation with respect to y. Recall that the Reynolds stress 
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interacts with the mean velocity gradient to increase or decrease the energy of pertur- 
bations. The distributions of Reynolds stresses for the cases under consideration are 
presented in Fig. 5. As the damping decreases Reynolds stress remains positive over a 
larger range showing increasing energy transfer from the mean flow 7 to the disturbances. 
Note also that the peaks are concentrated near the critical layer. 

The results presented in this paper demonstrate the applicability of the Chebyshev 
matrix collocation method to nonlinear eigenvalue problems in semi-infinite domains us- 
ing the companion matrix approach. The advantage of matrix collocation method in 
comparison with the tau method is the flexibility to use different coordinate transforma- 
tions and to impose different boundary conditions with relative ease. 
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FIGURE CAPTIONS 


1. Structure of C 0 after simple column operations. 

2. (a) Eigenfunction and (b) Derivative of eigenfunction for Re = 336 and u) — 0.1297; 
real (A) and imaginary (B) parts (Blausius Boundary Layer). 

3. (a) Eigenfunction and (b) Derivative of eigenfunction for Re = 598 and u> = 0.1201; 
real (A) and imaginary (B) parts (Blausius Boundary Layer). 

4. (a) Eigenfunction and (b) Derivative of eigenfunction for Re = 998 and u> = 0.1122; 
real (A) and imaginary (B) parts (Blausius Boundary Layer). 

5. Distribution of Reynolds stress; (A) Re = 336, u; = 0.1297 , (B) Re — 598, tu = 0.1201 
and (C) Re = 998, u; = 0.1122. 


Table 1: Comparison of the Eigenvalue Spectrum for Spatial Stability of Plane Poise uille 
Flow - Without Reduction (Re = 6000, w = 0.26) 



Bridges and Morris * 

Present Method 

Mode 

1 

1.00047-10.00086 

AT + 1 = 41 
1.00046-10.00086 

JV + 1 =51 
1.00047-10.00086 

2 

0.28323+10.02538 

0.28323+10.02538 

0.28323+10.02538 

3 

0.30165+10.04886 

0.30165+10.04886 

0.301 65+10.04886 

4 

0.31976+10.07532 

0.31976+10.07532 

0.31976+10.07532 

S 

0.33745+i0. 10492 

0.33748+10.10485 

0.33745+10.10492 

6 

0.35456+i0. 13782 

0.35664+10.13489 

0.35456+10.13782 

7 

0.37090+i0. 17425 


0.37089+10.17426 


Table 2: Comparison of the Least Damped Eigenvalue for Spatial Stability of Blausius 
Boundary Layer 


Re 

U3 

Jordinson 7 

Present Method 

336 

0.1297 

0.3084+10.0079 

0.30864+10.00799 

598 

0.1201 

0.3079-10.0019 

0.30801-10.00184 

998 

0.1122 

0.3086-10.0057 

0.30870-10.00564 
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